Bayesian estimation of pulsar parameters from gravitational wave data 
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We present a method of searching for, and parameterizing, signals from known radio pulsars in 
data from interferometric gravitational wave detectors. This method has been applied to data from 
the LIGO and GEO 600 detectors to set upper limits on the gravitational wave emission from several 
radio pulsars. Here we discuss the nature of the signal and the performance of the technique on 
simulated data. We show how to perform a coherent multiple detector analysis and give some insight 
in the covariance between the signal parameters. 

PACS numbers: 04.80.Nn, 02.50.-r 



I. INTRODUCTION 

Several kilometer-scale gravitational wave interferom- 
eters are now under construction or actively collecting 
data with unprecedented sensitivity 0, ■ To make the 
best use of these data, sophisticated analysis methods 
have been developed to search for astrophysical signals 
that are doubtless buried in the noise. A promising class 
of sources of gravitational waves are rapidly rotating, and 
structurally asymmetric, neutron stars. Several mecha- 
nisms have been proposed that could support a varying 
quadrupolar mass distribution in these neutron stars, and 
the subsequent continuous emission of gravitational ra- 
diation HQ. 

In this paper we present a detailed end-to-end descrip- 
tion of an analysis technique which we have developed 
to infer the parameters of such sources using data from 
interferometric gravitational wave detectors [5| . Here we 
test the performance of the technique on simulated data. 
Whether a signal is clearly present or not, the method 
can be used to set upper limits on emission strength. 

This method was successfully applied first to GEO 600 
and LIGO data from their first science run (SI) to set 
upper limits on the strength of gravitational wave emis- 
sion from pulsar J1939+2134 6]. The search was modi- 
fied and expanded to 28 isolated pulsars using data from 
LIGO's second science run (S2) 7]. This paper inves- 
tigates the methods used in these two papers and sets 
performance benchmarks for the algorithm. 

Searches for periodic gravitational wave signals from 
neutron stars are conventionally classed as blind, directed 
or targeted. A search is blind if no source parameters 
(such as sky position or spin evolution) are known a pri- 
ori, so that the size of the parameter space to be explored 
is maximal. Blind searches represent a computationally 
demanding problem, and highly efficient analysis tech- 
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niques must be used H,0|- Even with these methods, 
a fully coherent search using many months of data is 
computationally intractable, and the size of the param- 
eter space has the effect of decreasing the sensitivity of 
the search, as there is a good chance the noise will im- 
itate a relatively strong source somewhere in the space 
[Tfij . Directed (known location) and targeted (known lo- 
cation and phase evolution) searches have smaller num- 
bers of unknown parameters, and vastly smaller parame- 
ter spaces, making the detection level lower and increas- 
ing their sensitivity to gravitational waves. 

Radio pulsars are particularly interesting class of tar- 
geted sources because i) we can monitor their rotation 
and make a very good guess at the gravitational wave- 
form they produce and ii) their locations are known to 
high precision. In addition, the data pertaining to each 
pulsar can be restricted to a very narrow spectral win- 
dow. Real interferometric data is usually contaminated 
by a large number of instrumental spectral lines, so the 
effects of these can be reduced significantly by analyzing 
only those narrow bands containing pulsar data. 

This paper is structured as follows: sectionllTldescribes 
the nature of the gravitational wave signal that we are 
expecting from pulsars. Section II I II describes how we 
filter and greatly reduce the size of the data set using a 
semi-standard heterodyne technique. In Section IIVI we 
present the Bayesian methodology used for this analysis 
and describe the two likelihood functions used in previous 
work. We demonstrate the performance of the algorithm 
on simulated data in Section[V] Section lvTl concludes the 
paper with a brief summary. 



II. NATURE OF GRAVITATIONAL WAVE 
SIGNAL 

Here we summarize the form of gravitational waves 
emitted from a rotating rigid triaxial pulsar, described 
in detail in We take the special case of a triaxial 
ellipsoid rotating about its principal axis and therefore 
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emitting all its gravitational radiation at twice the ro- 
tation frequency. A freely precessing neutron star, with 
its spin and angular momentum axes non-aligned, would 
also emit at its rotation frequency but is expected to be 
strongly damped ^1 ■ The regularity of signals from the 
majority of radio pulsars suggests that most of them are 
not precessing on a short timescales, if at all. 

The gravitational wave amplitude from a triaxial neu- 
tron star seen from Earth is 



167T 2 GW r 2 

ho = ' A e, 



(1) 



where r is the distance to the pulsar, I zz its moment of 
inertia about the rotation (principal) axis, f r the rotation 
frequency of the pulsar, and e its equatorial ellipticity, 
defined in terms of its principal moments of inertia as 
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The gravitational wave strain on the detector will be 
frequency modulated, due to the relative motion of the 
Earth and the pulsar, and amplitude modulated by the 
antenna pattern of the interferometer. Following [81 Hlj. 
we can describe this measured signal as 



h(t) =^F+{t; ip)ho(l + cos 2 i) cos $(t) 
+ F x (t; ip)h cos 6 sin <!>(£), 



(3) 



where F + and F x are the antenna responses to the + 
and x polarizations respectively, ip is the polarization 
angle of the signal (determined by the position angle of 
the spin axis, projected on the sky), t is the inclination 
of the pulsar with respect to the line-of-sight, and Q(t) 
is the phase of the gravitational wave signal |2ij • 

We choose to time the signal phase evolution $(T) 
with respect to the solar system barycenter (SSB) , which 
is an inertial reference frame, so that to third-order in 
barycentric time, T, 



$(T) =0o + 2tt[/ s (T - To) + ~/ s (T - To) 2 
+ i/ s (T-T ) 3 ], 



(4) 



where <po is the phase of the signal at a fiducial time 
T , f s is the frequency of the signal (= 2/ r ), / s is the 
first frequency derivative, and / s is the second frequency 
derivative, all at time To. The transformation between 
the barycentric time (T) and the topocentric time at the 
detector (t) is 

T =t + St 

t ~t~ ^Roomer H - ^Shapiro ~t~ ^Einstein ^Binary: (5) 

where Ap{ ocmor is the classical Roemer delay, Ag Q is the 
Shapiro delay due to the curvature of space-time near the 
Sun, Aeq is the Einstein delay due to gravitational red- 
shift and time dilation, and Aeinary contains corrections 



related to the pulsar's orbit, which is zero for isolated 
neutron stars; see ^| f° r more details on these terms. 
However for pulsars in binary systems this term should 
include all the classical and relativistic corrections for the 
shifts in the time-of-arrival of the signal due to the mo- 
tion of the source within the binary system. We will not 
consider binary pulsars in this analysis, but for more de- 
tails on pulsar timing of binary systems see |l4l 1 1 5| and 
references therein. 

The second term in Equation [5] the Roemer delay, is 
the largest component (up to ~ 8.5 min) and due to the 
motion of the Earth within the solar system. In terms of 
the Earth's motion it is 



^Roomer 



r d k (r d -k) 2 



2cd 



(0) 



where r d is the position of the detector with regard to 
the SSB, k is a unit vector in the direction of the neutron 
star, c is the speed of light and d is the distance from the 
detector to the pulsar. In order to calculate the Roe- 
mer delay we need accurate knowledge of the position of 
the Earth with regard to the SSB. For our barycentering 
software |24| we use the solar system e phe merides pub- 
lished by the Jet Propulsion Laboratory [l6( . The second 
term in Equation is the timing parallax. This takes ac- 
count of the curvature of the wavefronts emitted from the 
source and is only significant for the closest sources. 

The Shapiro delay Ash ap i r o is a relativistic correction 
for the curvature of space-time near the SSB. Since this 
curvature is not negligible there will be an extra time 
delay in the arrival of a signal. In principle this delay can 
be as large as 120 (is for signals passing near the edge of 
the Sun and therefore becomes important for the analysis 
of signals from millisecond pulsars over periods of ~ 1 yr. 
The maximum contribution from Jupiter however is only 
200 ns and would not affect the sensitivity of a search. 

The Einstein delay describes the combined effect of 
gravitational redshift and time dilation due to the mo- 
tion of the Earth. This correction takes into account the 
varying gravitational potential experienced by a clock on 
the Earth as it follows its elliptical orbit around the Sun. 
Again this does not significantly affect signal searches. 



III. THE COMPLEX HETERODYNE METHOD 

Current ground-based interferometric detectors have 
broadband sensitivity to gravitational waves from fre- 
quencies of several tens of hertz up to several kilohertz. 
As a search for a continuous wave signal involves inte- 
grating for months or even years, the datasets involved 
can become very large indeed. However, the signal we 
are trying to extract in a targeted search is actually con- 
tained in only a very narrow frequency band, so accurate 
knowledge of the spin parameters of the source (from ra- 
dio or X-ray observations) allows us to reduce the size 
of this data set considerably. To do this, we perform a 
complex heterodyne, followed by filtering and resampling 
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of the data, to reduce its size by a factor of about ~ 10 6 
without loss of relevant information. Similar techniques 
have been applied in a wide range of optical, radio and 
gravitational searches for sinusoidal signals, most notably 
(in this context) |l7|. 

We choose to perform a complex, slowly evolving, het- 
erodyne on a targeted source to precisely unwind the ap- 
parent phase evolution of the source. However, for a sig- 
nal from a pulsar recorded by an interferometric detector, 
there is still a time varying component remaining in the 
heterodyned signal from the antenna response pattern of 
the interferometer. 

Since the source moves through the antenna pattern 
on a timescale that is much longer than the original peri- 
odic signal, after heterodyning we can re-sample the data 
with a much reduced sampling rate. In practice the new 
sampling rate is determined by our wish to monitor vari- 
ations in the interferometer noise floor, which changes 
on timescales of minutes, hours, and days. Since we keep 
both the real and imaginary part of the sample for each 
minute, our effective bandwidth is 1/60 Hz centered on 
the heterodyning frequency which is the instantaneous 
frequency of the signal at the detector. 

We take the calibrated output of a gravitational wave 
detector to be 

s(t) = h(t) + n{t), (7) 

where h(t) is a gravitational wave signal and n(t) is noise 
that is stationary over some time period but not neces- 
sarily Gaussian. Using Equation [3] we recast the signal, 
h(t), to 

h(t) = A{t)e 1 ^ + A* (i)e-**W , (8) 

where 

A(t) = ^F+(t;^)h {l + cos 2 i) - ~F x (t;ij)h cos t, (9) 

and A* is the complex conjugate of A. For a targeted 
pulsar we assume the frequency and frequency derivative 
terms are known from electromagnetic observations so 
that 

<f>{t) = 4(i + St) - 4>o (10) 

can be calculated to high precision. The heterodyning 
step involves multiplying the data from the interferome- 
ter by e -1 ^*-' to give 

Shct (t) = 8 (t)e-^ 

=A(t)e**° +>4*(t)e- < * , - 2i *W + n (t)e -< *W. 

(11) 

The heterodyning process removes the rotational phase 
evolution from the term in A, although this term will still 
vary over the day as the source moves through the an- 
tenna pattern of the interferometer. The second (upper 
sideband) term in A* will oscillate at nearly twice the 
gravitational wave signal frequency. 



We then apply a low-pass anti-aliasing filter to the het- 
erodyned data stream prior to averaging. We note that, 
because the original time series was real, we lose no inde- 
pendent information by rejecting the upper sideband. We 
use a series of three third-order Infinite Impulse Response 
(IIR) Butterworth filters to do this, with frequency cut- 
offs that can be adjusted to the characteristics of the 
data. The main requirement is to prevent spectral dis- 
turbances from outside our final 1/60 Hz data band being 
aliased into our calculation of the averaged data. The 
selection of the IIR filters and the sampling rate ulti- 
mately depends on the opposing needs to over-resolve 
the timescales on which the noise is non-stationary and 
for a narrow band to avoid nearby spectral lines. 

Finally we re-sample the filtered data to the post- 
filtering Nyquist rate and average the results (now s' het ) 
over a minute to form 

1 M 

where k is the minute index and M is the number of 
Nyquist samples in 1 minute (typically ~ 100). 

In practice there are computational advantages to per- 
forming the heterodyning, filtering and re-sampling pro- 
cess described above in two steps, starting with a fixed 
heterodyning frequency and a filter that reduce the sam- 
pling rate to about 4 Hz. A second (variable) hetero- 
dyne can then be performed on the data to remove the 
Doppler shifts due to the motion of the Earth. The ad- 
vantage is that the delay corrections between topocentric 
and barycentric time-of-arrival need only be calculated 4 
times, rather than (for LIGO and GEO) 16 384 times, 
per second. 

With the high frequency term in Equation El sup- 
pressed we have 

B k =^F+(t k ;i;)h {l + cos 2 tje** 

- *-F x {t k -^)h cos te** + n(t k )', (13) 

where n(tk)' is the heterodyned and averaged complex 
noise in bin k. By the central limit theorem, we would 
expect the noise, n(tfc)', to be well described by a Gaus- 
sian distribution, although the width of this distribution 
may change over time as the detector sensitivity evolves. 

The heterodyned gravitational wave signal in this re- 
duced data set depends on the same four unknown pa- 
rameters in Equations |21 and 0] ho, ip, (f>o, and t, which 
can be conveniently held as a vector, a. We proceed in 
the next section by calculating the (Bayesian) probability 
of the data given these parameters and finally, through 
the application of Bayes' Theorem and marginalization, 
we obtain posterior probabilities for each of these param- 
eters given the data collected. 
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IV. BAYESIAN FORMALISM 

We take a straightforward Bayesian approach for the 
following analysis and calculate the posterior probability, 
p{a\{B k \, J), of the pulsar parameters a given the binned 
data, {B k }. Bayes' Theorem tells us that 



p(a\{B k }J) = 



p(a.\I)p({B k }\8L,I) 
P({Bk}\I) 



(14) 



where a is the set of parameters inferred from data {B k }, 
given our model /, and with likelihood p({B k }\a, I). I 
remains constant throughout the analysis, and will be 
dropped from the following expressions to avoid clutter. 
It should of course be noted that all the inferences we 
make from the data assume this model to be true. Note 
that the posterior probability distribution given here as- 
sumes that a signal is present in the data. 

Our prior beliefs in the value of the parameters are 
held in p(a) and we will use the least informative priors 
for most of the parameters over their respective ranges. 
A change in polarization angle of 7r/2 on the sky is equiv- 
alent to a change of signal sign (i.e., a change in signal 
phase, 0o, of 7r), so consistent priors are </>o uniform over 
[0, 27r] and ?/> uniform over [— 7r/4, 7t/4]. The prior prob- 
ability density function for cos t is taken as uniform over 
[—1,1], corresponding to a uniform probability per unit 
solid angle for the orientation of the spin axis. 

The prior probability for ho is more interesting. In 
principle the prior for ho should reflect all our initial be- 
liefs on the gravitational wave strength, ho- If, for exam- 
ple we truly believe that gravitational wave emission is 
powered by the loss of kinetic energy from the pulsar (of 
known spindown rate), and that the moment of inertia of 
the pulsar is reasonably well constrained, then we should 
construct a prior that falls away sharply at strain lev- 
els which are above those consistent with this spindown 
upper limit. Of course we know that current detector 
sensitivities are insufficient to detect such a signal, and 
as a result the prior would overwhelm the broader like- 
lihood function. We would learn nothing new from the 
experiment since the posterior probability distribution 
function (pdf) would largely resemble the prior pdf we 
chose. 

At this stage in gravitational astronomy, a more useful 
statement would be concerned with what the observa- 
tions told us that was independent of spindown argu- 
ments, and therefore the prior should reflect this greater 
sense of ignorance. We cannot exclude the prior possi- 
bility of ho = 0, so a fully scale-invariant Jeffreys prior 
(oc 1/ho) would not be appropriate. However, we are 
interested in being able to set conservative upper limits 
on the strength of any signals, and this argues in favour 
of using a uniform prior for ho- A uniform prior favours 
larger values of ho over smaller values (e.g., the prior 
probability for the range 0.1 to 1 is 10 times less than for 
1 to 10), and represents, for most, an acceptable state 
of optimistic ignorance. The resulting upper limit for 
ho will therefore reflect the maximum value that could 



reasonably be though of as consistent with the data and 
has some additional merit because of that. In addition, 
a posterior based on a uniform prior for ho can be inter- 
preted as a (marginal) likelihood for ho and more easily 
incorporated into future analyses with other data. 

For further discussion choosing priors in cases, similar 
to this one, when the level of any signal may be below 
the sensitivity of the experiment, see |18j . Ultimately, if 
there is a strong detection the choice of the prior should 
not play an important role in the results since the like- 
lihood function would be sufficiently strongly peaked to 
define the posterior. Conversely, if no signal is present at 
the sensitivity level of the instrument, the prior takes on 
a greater significance. 

The full 4-dimensional posterior pdf contains all the 
information from our analysis but is difficult to interpret 
directly. It is therefore useful to reduce the dimensional- 
ity by marginalizing (integrating) over the less interest- 
ing (nuisance) parameters. The marginal distribution for 
one parameter can be viewed as a weighted average of all 
the distributions of this parameter given all the possible 
combinations of the other parameters. The parameter in 
which we are most interested is the gravitational wave 
amplitude ho, with a marginal pdf of 

p(ho\{B k }) ex JJJ p({B fe }|a)p(a)d0od^dcos t , (15) 

where the integral is performed numerically, over the full 
ranges of the nuisance parameters. The pdf for ho can 
then be normalised trivially. 

Even without a detection, placing upper limits on ho 
can be physically interesting as we are essentially con- 
straining the equatorial cllipticity of the neutron star. 
We define the upper limit of ho bounding 95% of the cu- 
mulative probability (from ho = 0) as the value /195 that 
satisfies 



0.95= / p{h \{B k })dh . 

lh =Q 



(16) 



Note that such a limit can be placed on ho even if most of 
the probability is to be found some distance from ho = 
and a strong signal is detected. 

In order to calculate the likelihood function we need 
to have a model of the signal in the processed data. The 
model of the signal that we are searching for in the data 
set, {-Bfc}, is obtained by processing the original gravita- 
tional wave signal h(t) in the same way that we processed 
the data to give 



y(t k ; a) =^F+(t k ; ^)h {l + cos 2 t )e**° 
- l -F x (t k ;i;)ho cosie** . 



(17) 



We note that the model is complex and that the only 
time varying component is the antenna pattern of the 
interferometers. The Nyquist frequency for this signal 
is well below our reduced sampling rate of one B k per 
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minute. In the following two sections we will present 
two different ways of evaluating the likelihood function 
depending on whether the variance of the data is known 
or unknown. 



A. Gaussian model — known variance 

Here we give the expression for the likelihood function 
assuming that we know, or can estimate accurately, the 
variance of the noise. We assume that the data comprises 
n samples of a signal, y(t k ) (see Eauation lT7|) . embedded 
in complex Gaussian noise N(0, a k ) of zero mean and 
known overall variance a\, so that 



B k =y k +N(0,a k ). 



(18) 



If the set of {B k } are independent, the likelihood of the 
data is simply the product of n bivariate normal distribu- 
tions. Note that the distribution is bivariate because the 
data are now complex. The real and imaginary parts of 
the Sfe's have independent noise components, each with a 
variance a\j2. The likelihood of the parameters is there- 
fore 



p({B fc }|a,{<7 fc }) - (V27ra fe )^ n exp - 



£ 

k=l 



\B k -y k \ 



(19) 

This Gaussian model for the likelihood was used for the 
first GEO 600 and LIGO analysis for signals from pulsar 
J1939+2134 6J. For this analysis, the noise level a k was 
not known but was estimated for each B k from the noise 
floor in a 4 Hz band of data around the signal frequency, 
assumed stationary for at least one minute. This gives 
240 x 2 points contributing to the estimate of the vari- 
ance, making the uncertainty in the point estimate of a k 
small enough to be ignored. The procedure is valid as 
long as the mean noise floor in the band is representative 
of the noise floor at the signal frequency, if there are no 
strong contaminating signals in the band and if the noise 
is sufficiently stationary. Although these requirements 
were largely met for 6], just one millisecond pulsar was 
involved in the study and they cannot be expected to be 
met in general. To address this an alternative model was 
developed for the S2 analysis Q. 



B. Gaussian model — unknown variance 

In the previous section we evaluated the likelihood 
function given the noise level, <j k , for each B k . Generally 
however, the noise level may not be known in advance 
or may not be well-estimated from the data. Here we 
described the likelihood function appropriate to this sit- 
uation, which was used in the analysis of the LIGO S2 
data 0. 

If a k is estimated from a tighter bandwidth, or over 
a shorter period, fewer data contribute and the uncer- 



tainty in its value may be too large to use a point esti- 
mate alone. Within our Bayesian framework the stan- 
dard (and correct) approach is to treat the noise level 
as another nuisance parameter and marginalise over it, 
without computing a point estimate at all |l9j . 

We begin by calculating the likelihood of a subset of rrij 
consecutive data points from {B k } which have a constant 
noise level aj. Once we have that expression, we will 
calculate the global likelihood simply using the product 
rule assuming that each segment of data is independent. 
We will again define n to be the total number of data 
points B k and let M be the number of segments of data 
that we have assumed have the same noise level, so that 



M 



(20) 



We can write the likelihood of the parameters, based on 
the jth subset of data and marginalised over aj as 

/■oc 

p({B k }j\a) oc / ^({Bfc}^, cr,- |a)dcr j 
Jo 

/•OC 

oc / p(a j \<i)p({B k } j \a ) a j )da j , (21) 
Jo 

where p((Tj |a) is the prior for the noise floor and the like- 
lihood p(aj\a)p({B k }j\a, aj) is given by Equation El As 
o j is a non-zero scale parameter we take a Jeffreys prior, 
uniform in log ay. 



p(cr,-|a) oc — (pj > 0). 

(To 



(22) 



Our final conclusions would be essentially unchanged if a 
uniform prior was used instead of this Jeffreys prior |20| . 
Here we assume that the <7j associated with each subset 
{B k }j is constant over the rrij samples. In other words, 
we assume that the noise level of the interferometer, in 
a narrow frequency band around the gravitational wave 
signal, is stationary for this time. However, we allow the 
noise floor to change between each subset of data {B k }j. 
This allows us to dynamically track the noise floor seen 
in real interferometric data, which will inevitably vary on 
some timescale, as the instrumental performance varies. 
The length of the jth subset, over the which the data is 
assumed stationary, can also be adjusted to reflect the 
known timescale of these variations. 

Using Equations ^| |^ and |22 the likelihood of the 
parameters based on a subset {B k }j of constant noise aj 
is 



P({Bfc»)cx / -n^-+T ex P - £ 



1 daj 

2(7? 1 J 



where rrij — k 2 (j) — feiy) + 1- This reduces to 



>({B fe }i|a) ck ( £ \B k ~y k f 

, k—k-, r.,-1 



(23) 



(24) 
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which is equivalent to a Student's i-distribution with 
2rrij — 1 degrees of freedom. Recall that the likelihood 
derived in Equation [2] is for a set of rrij data points 
Bk with the same <7j. The joint likelihood of all the M 
stretches of data, taken as independent, is therefore 



500 - 



400 - 



M I k 2(j) 

p({B fc }|a) cx J] ]T \B k -y k \ 

j \k=k 1(j - ) 



(25) 



j00 



We note that there is (of course) no explicit reference to 
the noise level, but the lengths of the stationary intervals, 
in j , can be adjusted to reflect the performance of the 
detectors. 



V. PERFORMANCE ON SIMULATED DATA 



A. Expected sensitivity 



J00 



100 - 








The analysis described above is optimal (in a Bayesian 
sense) for the data available from one or more science 
runs. It is however instructive to examine the long-term 
performance of the method on a large series of simulated 
data sets, both to confirm that the average performance 
complies with our expectations and to ease comparisons 
with methods based on sampling theory. 

To do this, we calculated the /195 upper limits from 
4000 simulated data sets, of length 10 days, varying the 
location of the putative source in the sky and the location 
of the detector in each set. The locations of the sources 
were picked randomly from a uniform distribution over 
the sky, and the detector locations were the GEO 600 and 
the two LIGO sites. 

From these we can express the average 95% upper limit 
(hgs} as a function of observation time T and single-sided 
noise power spectral density, S n (f). Empirically, from 
these simulations, 

(has) = (10.8 ± 0.2)y/S n (f)/T, (26) 

where the range accounts for the location of the detector. 
Figure^shows the distribution of hg 5 that contributed to 
this, for S n (f)/T = 1. Note that the width and skew of 
the distribution are relatively large, so the actual upper 
limit from an observing run could reasonably be up to a 
factor of two larger than (/195). 

B. Combining data from a network of detectors 

Several gravitational wave detectors are currently col- 
lecting data, and ideally we should be able to use the 
observations from all detectors in a coherent manner in 
order to draw the best possible inference about the source 
parameters. In a Bayesian analysis all observations enter 
via the likelihood function. Assuming that the noise from 
each interferometer is independent, by the product rule 



FIG. 1: Distribution of 95% upper limits on ho for 4 000 
simulations, using sources randomly located on the sky with 

Sn(f)/T = 1. 

the global likelihood is simply the product of the individ- 
ual likelihoods. For example, by combining observations 
from GEO 600 and the three LIGO interferometers, we 
would get 

4 

p({B k } Joint I a )=Hp({Bk}M (27) 

8=1 

where the product is over the 4 km Hanford interferome- 
ter (HI), the 2km Hanford interferometer (H2), the 4km 
Livingston interferometer (LI), and GEO 600. 

This likelihood contains all the information on the 
source parameters that is contained in the data, opti- 
mally combining the data from all the interferometers in 
a coherent way. Note that the observation periods can be 
different and so can the sensitivity of the detectors, al- 
though for detectors with very different sensitivities, this 
will closely approximate the likelihood based on just the 
most sensitive instrument. 

For illustration, we generated four sets of data of 10 
days length with Gaussian noise (fi = and a = 1) as 
if from GEO 600 and the three LIGO interferometers. 
As these IFOs are modelled as having the same sensi- 
tivity, we would expect the coherent results to be ap- 
proximately s/Z times tighter than the individual results 
(distinguished from a factor of four increase in observ- 
ing time only by the differing antenna patterns of the 
instruments). The four posterior pdfs for each detector 
as well as the joint multi-detector posterior pdf for ho 
are shown in Figure [3 The individual 95% upper limits 
are 0.15 for GEO 600, 0.16 for HI, 0.18 for H2, and 0.13 
for LI giving an average of 0.155. The combined 95% 
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FIG. 2: Multi-detector posterior pdfs with simulated data. 
The solid line represents the joint marginalised posterior pdf 
for ho using the data from four separate interferometers. The 
dashed lines are the corresponding pdfs from the individual 
detectors. 



upper limit, on the other hand, is 0.08, which is indeed 
approximately a factor of 2 better than the average of the 
limits from the individual detectors. This technique was 
first applied to real gravitational wave data in Ml. The 
equivalent multi-detector analysis using the jF-statistic 
method |Sj has recently been developed by Cutler and 
Schutz [2lf . 

It is important to realise that the posterior curve de- 
rived from a particular observation represents a prob- 
abilistic statement about the value of ho based on the 
data in hand and may, if we are unlucky, be wildly at 
odds with the truth. As a result the upper limit derived 
from one instrument alone will occasionally be lower than 
that from the coherent combination of instruments. 



C. Effects of changing the noise level 

The widths of the marginal posteriors depend on both 
the level of the noise and the covariance of the param- 
eters. Here we demonstrate the noise dependence by 
analysing three sets of data with different, sometimes 
modulated, noise variances. Each data set corresponds 
to 10 days of observations. The first contains Gaussian 
noise with /i = and a = 1. For the second data set, the 
noise level alternates each 30 minutes between a = 10 
and a = -^j^- For the third data set, the noise level al- 
ternates each 30 minutes between a = 100 and a = —K=. 

5v 2 

Two time series plots showing representative stretches of 
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FIG. 3: Time series of Bk showing the effect of a changing 
noise level. The dots are the \Bk\ and the line represents the 
variance of the -Bfc's. The top figure represents the first data 
set with constant variance and the bottom figure represents 
the second data set with alternating variance. 

data from the first and second sets are shown Figure 

For this test, we repeated and averaged the posterior 
pdfs for 100 generations of the data sets described above. 
The average marginalised posterior pdfs for ho are shown 
in Figure Using the 66% upper limit on ho to charac- 
terise the width of the pdfs we have ho < 0.095 for case 
1, h < 0.047 for case 2, and h < 0.019 for case 3. The 
second posterior is narrower than the first by a factor 
of 2.02 and the third by a factor of 5.05. For no signal, 
these results are roughly what we would expect: in the 
two cases with alternating noise levels, about half of the 
data has very low sensitivity compared to the other half, 
and we can assume this noisy half does not play a sig- 
nificant part in the posterior. Therefore these cases are 
approximately equivalent to a continuous observation at 
the greater sensitivity level but for half the full observ- 
ing period, and the reduced time lowers the sensitivity by 
•s/2- Thus compared to the first case with constant noise, 
we would expect the two sets with alternating noise levels 
to have widths which are narrower by factors of about 2 
and 5, which is what we see. 



D. Covariance between parameters 

To illustrate the covariance between the signal parame- 
ters we generated a 10 days time series containing a signal 
with the following parameters and Gaussian noise of unit 
variance: h = 0.25, cost = 0.1, (j) = 180°, ip = 0.0°. It 
is clear from the emission model that if cos i ^ 0, ho and 
cost are strongly anti-correlated, as are ip and <po- The 
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FIG. 4: Averaged marginalised posterior pdfs for ho for three 
scenarios: constant unit variance (solid line), alternating noise 
level between a — 10 and a — (dashed line), and alter- 
nating noise level between a — 100 and a = (dotted 
line). 

correlation can be seen clearly in the probability den- 
sity contour plots in Figure [SJ The covariance between 
cos i and ho contributes strongly to the overall width of 
the marginal pdf of ho, making the precise value of ho 
somewhat difficult to determine even under conditions of 
relatively high signal-to-noise ratio. The correlation be- 
tween ho and cos t is also visible when no gravitational 
wave signal is injected in the data (ho = 0), largely be- 
cause of our choice of a uniform prior for ho which holds 
out high hopes for a strong signal in the data (Figure EJ. 
When no such signal is seen, this is interpreted as an in- 
dication that the pulsar is oriented unfavorably and the 
posterior probability slightly increases around cos t = 1, 
where only the '+' polarization is present. 



VI. CONCLUSIONS 

In this paper we have presented an end-to-end 
Bayesian method of searching for, and parameterizing, 
gravitational waves from known pulsars. The method 
involves processing the raw data to reduce the number 
of samples required for the analysis. We calculated the 
likelihood function for given model parameters from the 
decimated data, so reducing computational requirements. 
The algorithm has been validated by retrieving the cor- 
rect signal parameters from simulated data. We have also 
shown than it is easily adapted to deal with a network of 
detectors. 

This methodology was initially developed for targeted 
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FIG. 5: Equally-spaced contours of constant probability den- 
sity for (top) the joint posterior probability of cos t and ho and 
(bottom) ip and <f>o. The marker indicates the location of the 
simulated signal. 

searches with known locations and spin evolutions of 
the sources. Further work has now been done study- 
ing the feasibility of expanding the numbers of parame- 
ters by taking advantage of Monte Carlo Markov Chain 
(MCMC) techniques [22|. These techniques are required 
when the number of unknown parameters is significantly 
increased, when the method presented here would be too 
computationally intensive. In a future paper we will ad- 
dress how this method has been adapted to search for 
gravitational emission from pulsars in binary systems. 
The algorithm presented in this paper, with the binary 
modification, is currently being applied to GEO 600 and 
LIGO data from the S3 and S4 science runs. 
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FIG. 6: Equally-spaced contours of constant probability, 
showing the covariance between ho and cost, for Gaussian 
noise with no signal but a uniform prior for ho. 
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